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Shear flow significantly affects the transport of swimming algae in suspension. For 
example, viscous and gravitational torques bias bottom-heavy cells to swim to- 
wards regions of downwelling fluid (gyrotaxis). It is necessary to understand how 
such biases affect algal dispersion in natural and industrial flows, especially in view 
of growing interest in algal photobioreactors. Motivated by this, we here study 
the dispersion of gyrotactic algae in laminar and turbulent channel flows using di- 
rect numerical simulation (DNS) and the analytical swimming dispersion theory 
of Bees & Croze [l[ . Time-resolved dispersion measures are evaluated as functions 
of the Peclet and Reynolds numbers in upwelling and downwelling flows. For lam- 
inar flows, DNS results are compared with theory using competing descriptions 
of biased swimming cells in shear flow. Excellent agreement is found for predic- 
tions that employ generalized- Taylor-dispersion. The results highlight peculiarities 
of gyrotactic swimmer dispersion relative to passive tracers. In laminar downwelling 
flow the cell distribution drifts in excess of the mean flow, increasing in magnitude 
with Peclet number. The cell effective axial diffusivity increases and decreases with 
Peclet number (for tracers it merely increases). In turbulent flows, gyrotactic effects 
are weaker, but discernable and manifested as non-zero drift. These results should 
significantly impact photobioreactor design. 

Keywords: algae, swimming microorganisms, Taylor dispersion, DNS, 
turbulent transport, bioreactors 

1. Introduction 

In natural bodies of water and in industrial bioreactors, microscopic algae expe- 
rience laminar and turbulent flows that play a critical role in their dispersion, 
proliferation and productivity @; Q • The dispersion of passive tracers in fluid flows 
is well understood, particularly in the industrially relevant cases of flows in pipes 
and channels G. I. Taylor first realised that the dispersion of passive tracers in 
a pipe, caused by the combination of fluid shear and molecular diffusion, could be 
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described in terms of an effective axial diffusivity [B[ . A similar result also holds for 
turbulent pipe [6f] and channel Q flows. Taylor's pioneering analyses of dispersion 
in a pipe inspired a series of studies that placed the understanding of the disper- 
sion of neutrally buoyant tracers on a firm footing 0; Q . Particles whose density 
differs from the suspending medium exhibit more complex behaviour (for example, 
they can accumulate at walls in turbulent channel flows; see 1^; U| and references 
therein) . 

Swimming single-celled algae are known to respond non-trivially to flows. For 
example, the mean swimming direction of biflagellates, such as Chlamydomonas and 
Dunaliella spp., is biased by imposed flows 12; T^. This bias, known as gyrotaxis, 
results from the combination of viscous torques on the cell body, due to flow gradi- 
ents, and gravitational torques, arising from bottom-heaviness and sedimentation 



14| . In the absence of flow gradients, the gravitational torque leads cells to swim 



upwards on average (gravitaxis) . Recent simulations and laboratory experiments 
have shown how inertia and gyrotaxis can play significant roles in the formation 
of patchiness and/or thin layers of toxic algae in the ocean [2 E3|- Gravitaxis 
and gyrotaxis can also couple in a complex fashion to other biases due to chemical 
gradients (chemotaxis) and light (phototaxis (HI). Therefore, the fact that cells 
can actively swim across streamlines and accumulate in specific regions of a flow 
forewarns that the resultant dispersion of biased swimming algae in a complex flow 
field is nontrivial. 

In a recent study, Bees & Croze extend the classical Taylor-Aris analysis of 
dispersion in a laminar flow in a tube to the case of biased swimming algae [l[. 
The Bees & Croze dispersion theory provides a general theoretical framework to 
describe the dispersion of biased micro-swimmers in confined geometries, such as 
pipes and channels. To make predictions for the dispersion of particular organisms, 
specific functional expressions for the swimming characteristics need to be obtained 
from microscopic models. For example, using expressions for the mean swimming 
orientation q and diffusivity D from the so-called Fokker-Planck model |13| , Bees 
& Croze found that gyrotaxis can significantly modify the axial dispersion of swim- 
ming algae (such as Chlamydomonas spp.) in a pipe flow relative to the case of 
passive tracers (or dissolved chemicals) More recently, Bearon et al. [HI ob- 
tained predictions for swimming dispersion in laminar flows in tubes of circular 
cross-section by using an alternative microscopic model known as generalised Tay- 
lor dispersion (GTD), formulated for swimming algae 2(J 21 1 and bacteria 22 1. The 
GTD model is considered superior, as it incorporates correlations in cell positions 
due to cell locomotion within local a flow, as well as the gyrotactic orientational 
biases considered in the Fokker-Planck description. The qualitatively distinct pre- 
dictions, however, have yet to be tested. Direct tests of the models might involve 
microscopic tracking of isolated cells swimming in prescribed flows. The Bees & 
Croze theory allows the predictions of the two models to be tested for entire popu- 
lations in macroscopic experiments [23j j and individual-based numerical simulations, 
as we shall demonstrate in this paper. 

Microalgae currently play a prominent role in biotechnology. For example, they 
are cultured commercially as nutritious fish and crustacean feed in aquaculturc, 
and for high- value products that they can synthesize, such as /3— carotene j24[. 
Microalgae also hold significant potential for future exploitation; they can pro- 
vide a sustainable biofuel feedstock that does not compete for arable land with 
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food crops. Despite significant efforts, however, current production systems remain 
too inefficient for algal biofuels to be commercially viable (25|. Algal culture com- 
monly occurs in open or closed bioreactors, with production systems employing 
three stages: cell culture, harvest and downstream processing. Open bioreactors 
consist of one or more artificial ponds, quiescent or stirred, where stirring allows 
to obtain greater concentrations per reactor area [24 |. A common design is the 
raceway pond, in which rotating paddles generate flow (sec figure [TJi) . Typically, 
raceway ponds are shallow (depths L = I — 30 cm), which aids cell exposure to 
light. Typically, the flow is turbulent (speeds of 30 cm/s), facilitating light expo- 
sure and gas exchange. Pond bioreactors are relatively cheap to operate, but are 
limited by their low areal productivity and susceptibility to contamination. This 
restricts the range of viable species in open pond culture to robust, high- value ex- 
tremophiles, which grow in prohibitive conditions such as high salinity or low pH. 
A greater variety of fast-growing algae can be cultured in closed bioreactors, where 
single species are propagated within sealed transparent containers. Common ge- 
ometries are flat, cylindrical (columnar or tubular), and annular, see figure [T|>c. 
As for raceways, the constraint of light penetration limits the smallest dimension 
to fO cm. In tubular bioreactors the flow-speed typically is 50 cm/s, whilst flat 
panel and columnar /annular set-ups operate at 1-10 cm/s [2f|. Although closed 
bioreactors allow for rapid growth conditions with small footprints, current designs 
are relatively expensive to manufacture, operate and maintain, and the algae are 
susceptible to invasion by other less-useful species. 




Figure 1. Examples of photobioreactors: (a) raceway pond (top view), showing the pad- 
dlewheel and guiding baffles; (b) tubular array, where suspensions are typically driven by 
a pump; (c) air-lift flat panel (3D view) and draft-tube (cross-sectional side view), where 
rising bubbles generate flow as well as providing mixing and aeration. The regions high- 
lighted with broken lines in the top diagrams are reproduced below to indicate the length 
scale L and magnitude U of typical flows. For the airlift in (c), both upwelling U u and 
downwelling Ud flows are shown. 



Existing bioreactors operate under both laminar and turbulent flow conditions. 
The transition between the two flow regimes depends on the ratio of inertial and 
viscous forces in a fluid. This is expressed by the dimensionless Reynolds number 
Re = UL/i/, where U and L are the characteristic speed and length scale of the 
flow, respectively, and v is the suspension kinematic viscosity 27[ . In general, flows 



Article submitted to Royal Society 



4 



O. A. Croze et al. 



with Re > 2000 will be turbulent, though the transition to turbulence depends on 
geometry and particular flow conditions. Approximating the kinematic viscosity by 
that of water and using the flow speed and length scales above, we see that flows in 
air-lift reactors can be both laminar or turbulent (Re> 100), while raceway ponds 
and tubular reactors are always turbulent (Re> 30000). The mixing properties of 
turbulent flows are thought beneficial for cell growth, although clear evidence for 
this is lacking. With air-lifts the mixing is caused by rising bubbles, which also 
provide aeration, and the flow does not need to be turbulent. Both cells and nutri- 
ents will disperse in the above flows, with cell growth and productivity depending 
critically on such dispersion. 

Current photobioreactor designs assume that cells disperse like tracers, whether 
they swim or not j28[. However, as detailed above, recent studies show that the dis- 
tribution and consequent dispersion of swimming cells in flows, including biotech- 
nological interesting species like Dunaliella, should be very different from that of 
passive tracers [H Il9j . For example, if cells disperse differently to nutrients in a 
reactor they will separate, which could have catastrophic consequences for growth. 
Inspired by the possibility of taking advantage of the peculiarities of swimmer re- 
sponse to flow to improve bioreactor operation, we study here the dispersion of 
gyrotactic swimming algae in laminar and turbulent channel flows. We carry out 
time-resolved direct numerical simulations, comparing results to predictions from 
the Bees & Croze dispersion theory applied to channel geometry. The paper is 
structured as follows. In Section 2, we present the simulation methods and outline 
a derivation of the theory for the new geometry. In Section 3, we present simulation 
results for the dispersion of biased swimming algae in laminar and turbulent flows, 
comparing theory and simulation results for laminar flows. In Section 4, we inter- 
pret the results, and discuss their implications for dispersion in photobioreactors 
and the environment. 



2. Methods 

(a) Direct numerical simulations (DNS): governing equations 

We simulate the dynamics of a population of biased swimming micro-organisms, 
typically N = 2 x 10 5 and 10 6 individuals per simulation (see table [J), placed 
in laminar and turbulent flows. Each microswimmer is modelled as a spheroidal 
particle whose position Xj, i — 1..N, evolves according to 

f = u + V.p (2.1) 

where V s is the mean cell swimming speed, u is the local fluid velocity and p is the 
local particle orientation. 

The orientation p of each swimmer evolves in response to the biasing torques 
acting upon it. Making the realistic simplifying assumption that swimming cells 
have a spheroidal geometry, the reorientation rate of the organisms is defined by 
the inertia- free balance of gravitational and viscous torques 29|; 13 1 



^ = -L [k - (k ■ p)p] + x p + a (I pp) E p + IV, (2.2) 
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where k is a unit vector in the vertical direction; B = p,a±/(2hpg) is the gyrotactic 
reorientation time (h denotes the centre-of-mass offset relative to the centre-of- 
buoyancy, a±_ is the dimensionless resistance coefficient for rotation about an axis 
perpendicular to p, p and p are the fluid density and viscosity, respectively); ao = 
(a 2 — b 2 )/(a 2 + b 2 ) is the eccentricity of the spheroids with major axis a and minor 
axis b; and, finally, E is the rate of strain tensor and uj the vorticity. The noise T r 
is added to simulate the stochastic rotational diffusivity of a swimming cell. In the 
simulations it is implemented with random angular steps of magnitude ^2d r , so 
that cells diffuse with rotational diffusivity d r [3(J [3l[ . 

The flow field u is obtained by solving the Navier-Stokes equations for an in- 
compressible viscous fluid, such that 

— - = --\7p + v\7 2 u, V-u = 0, (2.3) 

Dt p 

where ^ = + (u • V) denotes the material derivative. Here, for simplicity, we 
shall assume cells to be neutrally buoyant, neglecting feedback from the particles 
to the flow (these effects are explored elsewhere [l|). 



(b) DNS: geometry, scaling and statistical measures of dispersion 

We consider a standard channel geometry: two flat plates parallel to each other, 
infinite in extent, and separated by a gap of size 2H. When focusing on swimming 
cells it would seem natural to rescale length by the plate half-width H and time 
by H 2 d r /V 2 , where V s is the mean swimming speed and d r is the cell rotational 
diffusivity, the characteristic time a cell swimming with V s takes to diffuse across 
H with diffusivity V 2 d~ Y . However, the parameter values obtained from this 'cell- 
based' rescaling are too large and thus numerically inconvenient for simulations. 
Thus we adopt a 'flow-based' rescaling in terms of the characteristic length H and 
the flow-based timescale H/U c , the time taken for a flow with centerline speed U c 
to advect a cell by H. In terms of this rescaling, the dimensionless equations of 
motion for a biased swimming cell are 



u*+v s p, (2.4) 



dt* 

= ^ 1 [k-(k.p)p] + ^* xp + ao (I-pp)-E*.p + r;, (2.5) 

= -V> + i?e" 1 V* 2 u*, V-u*=0, (2.6) 

where starred quantities denote dimensionless variables. In particular, we define the 
dimensionless swimming speed v s — V s /U Cl the gyrotactic parameter 77 = 2BU C /H 
and the Reynolds number 

Re=^, (2.7) 

where U c is the centerline speed, H is the channel half-width and v the fluid's kine- 
matic viscosity. Non-dimensionalising the noise T r also defines the dimensionless 
rotational diffusivity d* — d r H/U c . 

We adopt a Cartesian coordinate system where the mean flow in the x (stream- 
wise) direction, varies with the wall-normal coordinate y, and is independent of z 
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(spanwise direction). We integrate equations (12.41) . (I2.5[) and (I2.6[) numerically (see 
supplementary materials) to find x?(t*) = [x* (i*), y*(t*), z*(t*)}. Enumerating the 
number of cells N(x) at a given position x in a bin of fixed volume AV, we obtain 
the cell concentration c(x) = N(x)/AV. From the cell coordinates can we further 
define statistical measures of the cell dispersion in a flow: the adimensional drift Ag 
with respect to the mean flow U — (2/3)C/ c ; the effective streamwise diffusivity D* 
and the skewness of the cell distribution, 7. These measures of dispersion are given 
by 

. */ *s dm! 2 „. 1 1! „ , . „. m% — 3m!m% + 2m! 3 

(2.8) 

where m* = ^5Zj(x*) p (P = 1>2,3) are the distribution moments and Var(ai*) = 
m| — to* 2 is the variance of the cell distribution. The statistical measures can be 
transformed from the flow-based scaling with characteristic time-scale H/U c to the 
cell-based scaling with time-scale H 2 d r /V 2 by the transformations 

t* ->■ t = t*/Pe, A* ->■ A = A* Pe, D* D c = D* Pe, (2.9) 

where Pe is the cell Peclet number defined with respect to the center-line speed (see 
equation (|2.13[) ). Note that the skewness 7 does not depend on the scaling used. 
We shall present results in terms of the cell-based scaling and compare their limit 
at long times with analytical predictions. 



(c) Analytical theory: dispersion in laminar flows at long times 

Here, we shall obtain analytical predictions for the dispersion of algae swim- 
ming in laminar channel flow that will be compared with the results from the 
simulations. The general Bees & Croze [l[ continuum dispersion theory for biased 
swimmers will be applied to the new channel geometry; it is valid in the long-time 
limit (t ^> H 2 d r /V 2 ). The derivation for channel flow is similar to the pipe flow 
example in Bees & Croze, so we provide an outline here. Those readers that are 
less mathematically inclined may skip the derivation of the results in this section. 

We shall begin with the continuity equation for the cell number density n: 

On 

— = -V-[n(u + V c )-D-Vn]. (2.10) 

The flow velocity u is the solution of the Navier-Stokes equation (|2.3p . We adopt 
the same coordinate system described above for the DNS. For laminar flows such 
that the flow downstream is only a function of the wall-normal coordinate y, the 
flow velocity can be expressed as 

u{y)=u(y)e x = U[l + X (y)}e x , (2.11) 

where U is the mean flow speed and x(y) describes the variation in the streamwise 
direction about this mean. It is convenient to translate to a reference frame travel- 
ling with the mean flow, and non-dimensionalise lengthscales by H and timescales 
by the time to diffuse across the channel H 2 d r /V 2 . Thus, x = (x — Ut) /H, y = y/H, 
z = z/H and t = V 2 t/(H 2 d r ), where hats denote dimensionless variables. 
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We assume unidirectional coupling of the cell dynamics to the flow (for bidi- 
rectional coupling due to non- neutrally buoyant cells see cells are biased by 
shear in the flow, so that the cell swimming velocity V c = V s q, where q is the 
mean swimming direction, and diffusivity tensor D are functions of local flow gra- 
dients. Here, we consider analytical predictions for spherical cells (olq = 0), so cell 
orientation is only a function of vorticity, u> = V x u = — x' e z- 

Consider planar Poiseuille flow, \ — (1 — 3y 2 )/2 and U = (2/3) t/ c |32| . where 
U c is the centerline (maximum) flow speed. Equation (|2.10j) becomes 

Ufl 

— = V • (D • Vn) — (2/3)Pe X n x - /3V • (nq) , (2.12) 

with 

U c Hd r Hd r ,„.„> 

Pe=— — — , and (3 = — — , (2.13) 

where hats have been dropped for clarity. Pe is a Peclet number, the dimensionless 
ratio of the rates of transport by the flow and swimming diffusion, and j3 is the ratio 
of channel half-width to the length a cell swims before reorienting significantly. 
Alternatively, we can re-write f3 = HVgd r /V^ and interpret it as a 'swimming' 
Peclet number, the ratio of the rates of transport by swimming and diffusion. No- 
flow and no-flux boundary conditions will be applied to (12 . 12|) . such that 

u = and n ■ (D • Vn - f3qn) = 0, on E, (2.14) 

where n is the unit vector normal to the channel boundary E. 

As the flow is translationally invariant along x, the mean swimming direction 
and diffusivity tensor are independent of x: D — D(y), q = q(y). This permits a 
treatment of dispersion using moments in a similar vein to that of Aris Q . The pth 
moment with respect to the axial direction is defined as 

/+oo 
x p n(x,y,t)dx, (2-15) 
-oo 

provided x p n{x,y 1 t) — > as x — > ±oo. We denote cross-sectional averages by 
overbars. The cross-sectionally averaged axial moment is thus 



1 



m pW = c p = 2 J c p{y^) d v- ( 2 - 16 ) 

In Cartesian coordinates equation (|2.12[) becomes 

n t = iiy'n,, -U, :l n ■ /)"'//,;, • /;"'//,, (2.17) 
-[(2/3)Pex + (3<f]n x + D xx n xx . 

Multiplying by x p and integrating over the length of an infinite channel we obtain 
the moment evolution equation 

c P ,t = (D yy c p . y - {3qyc p - P D x yc p ^) y - P D xy c p _ hy (2.18) 
+p[(2/3)Pe X + (3q x }c p ^ + p(p - l)D xx c p - 2 . 
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Averaging over the cross-section and applying the no-flux boundary conditions 
(gUJ) , yields 



m Ptt = p(p - 1)-D*%_ 2 - pD^Cp-Ly + p[(2/3)Pex + Pq x ] c v . x (2.19) 

from which we calculate measures of cell dispersion. The drift above the mean flow, 
A = Hindoo ^TOi(t), is given by 



A = -D*vY °> + [(2/3)Pe X + /9?»] Y °, (2.20) 

where 



F °(y) = exp^£-|M y d.)^oxp( ;/ ^— ,/* ) ) (2.21, 



is the zeroth axial moment (normalised concentration profile) . Similarly the effective 
diffusivity, D e = limine [m 2 (t) - m\(t)\ , is given by: 



De = -Dwg' + [(2/3)Pex + Pq* - A ] g + D XX Y§, (2.22) 

where g(y) = Y° /» (^- kdg^ffi )^ with A (y) = [-ZW °'+[(2/3)Pe X + /3<f ] Y °]ds 

and Too (y) = J Iq ^ 5, The weighting function ^(y) controls the drift and g(y) 
(related to the first axial moment) controls the value of the diffusivity. 

To make predictions from (|2.20l) and (12.221) we require expressions for D and 
q from microscopic models of the statistical response of cells to flow. Two main 
models have been proposed for biased swimming algae: the Fokkcr-Planck (FP) 
model [Ullil and generalized Taylor dispersion (GTD) USUI El]. For s P herical 



cells, each model predicts how nondimensional swimming direction q and diffu- 
sivity D depend on two nondimensional quantities: a(y) — —x'(y)U/(2d r H) = 
— x'(y)(l/3)(Pe//3 2 ), the ratio of reorientation time by rotational diffusion to that 
by shear (vorticity), and A = l/(2d r B), the ratio between the time of reorientation 
by rotational diffusion (l/2)d~ 1 and the gyrotactic reorientation time B. The func- 
tional forms for the transport parameters q(er) and D m (cr), where subscripts m =F 
and G denote the FP and GTD models, respectively, are given in appendix A. In 
particular, the two models differ qualitatively in their predictions for the diffusiv- 
ity as a function of the shear rate and therefore they provide different predictions 
for cell distribution and dispersion. The dispersion predictions (12.20)) and (|2.22[) 
were evaluated numerically with Matlab (Mathworks, Natick, MA) using q(er) and 
D m (<r) from the two models, see supplementary materials. 



(d) Parameters, simulation time and flow profiles 

In simulations and theoretical evaluations we used the following cell parameters 
based on C. augustae: V s = 0.01 cm s _1 (swimming speed), d r = 0.067 s _1 (cell 
rotational diffusivity), and B — 3.4 s (gyrotactic reorientation time) [I]. Further, 
we consider the flow of suspensions taken to have the same viscosity as water, 
v = 0.01 cm 2 s _1 . The cell eccentricity ao is also held fixed. The analytical theory 
presented above assumes that cells are spherical, ao = 0, but laminar and turbulent 
simulations have also been performed for elongated cells with ao = 0.8. With these 
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Re (Lam/Turb) 


Pe 


P 


ao 


4(L) 


27 


1.34, 6.7, 33.5 





100 (L) 


670 


1.34, 6.7, 33.5 


(for all P), 0.8 (for /3 = 6.7) 


250 (L) 


1675 


3.35, 33.5 





2500 (T) 


16750 


33.5 





4200 (T) 


28140 


55 


0, 0.8 


10000 (T) 


67000 


67 






Table 1. Nondimensional parameters for the different cases considered in the DNS and in 
the analytical model predictions for the laminar case. Re is the Reynolds of the flow, Pe 
and ft are flow and swimming Peclet numbers, respectively, defined by equations i2.13\) , 
and qq is the cell eccentricity (> for elongated cells, zero for spheres). 



parameters fixed, choosing the centreline speed U c and channel width H gives the 
nondimensional flow-based parameters used in the DNS: v s = V s /U c (dimensionless 
swimming speed), r\ — 2BU C /H (gyrotactic parameter), d* = d r H/U c (dimension- 
less rotational diffusivity) , and Reynolds number Re= U c H/v. For the test run 
with passive tracers an additional noise term was added to equation [23] to simulate 
the translational diffusivity D t , nondimensionalised as D* t — D t / \U C H). These flow- 
based parameters can be transformed into cell-based nondimensional parameters for 
comparison with analytical predictions. From the definitions above and equation 
(|2.13[) it can be shown that Pc= d*/v^ (flow Peclet number), j3 = d*/v s (swimming 
Peclet number) and a(y) = — (l/3)x'(y)/d* (local dimensionless shear-rate). Since 
x'(y) — — 3y for plane Poiscuille flow, the maximum dimensionless shear is given 
by |<7maa;| = = TjX, where we recall A = l/(2Bd r ), the nondimensional bias 

parameter. Simulations are more readily interpreted and compared to analytical 
theory in terms of these parameters (shown in table [J) . 

The dimensionless flow profiles corresponding to the Reynolds numbers used in 
this study are plotted in figure[2]for the benefit of the reader. Laminar flows are self- 
similar, so the nondimensional flow has the same parabolic profile for all Re. Time- 
averged turbulent flows have distinct profiles that depend on the Reynolds number. 
Note that the number of degrees of freedom in the turbulent flow simulations scales 
as Re 9 27| . This makes large Re simulations computationally expensive. For this 
reason we do not investigate dispersion for flows beyond Rc= 10000. 



3. Results 

(a) Passive tracers: classical dispersion 

The dispersion of passive tracers, such as molecular dyes or nonmotile cells, is 
generally well understood. In laminar channel flow passive tracers are transported 
on average at the mean flow speed; there is no drift relative to mean flow: Ao = 0. 
The effective axial diffusivity D e , is given at long times by the Taylor- Aris result 
D e = l+ (8/945)Pe 2 [j^H- As a benchmark test, we carried out direct numerical 
simulations for passive tracers (solving equation (|2.4p with v s — and translational 
noise to simulate molecular diffusivity, see methods). A typical result is shown in 
figure [31 for Pe= 3000. As expected at long times Ao — > and the Taylor- Aris 
diffusivity prediction, D e — 76191, compares very well with simulation results. The 
skewness slowly tends to zero at long times as t~ 5 (see supplementary material), 
suggesting approach to a normal distribution 0. The above results for passive 
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Figure 2. Flow profiles for the dimensionless flow speed as a function of y for the laminar 
cases (self-similar for all Re) and time-averaged flow profiles for the turbulent cases (Re, 
Pe)= (2500, 16750), (4200,28140) and (10000,67000), as shown. 




Figure 3. Time dependence of the effective diffusivity, D e , for passive tracers in a laminar 
flow with Pe= 3000. For long times D e ~ 76191, the constant value predicted from the 
classical Taylor- Aris dispersion for passive tracers (see text). The inset shows the expected 
zero drift, Ao, above the mean flow and the skewness, 7, tending to zero at long times, 
suggesting an approach to the expected normally distributed axial concentration profile 



tracers will be compared with the dispersion of gyrotactic swimmers. They can also 
be used to check our analytical theory in the limiting case of unbiased swimmers 
with no gyrotactic response to flow. In this limit, q = and the diffusivity tensor 
is isotropic D = 1/6, where I is the identity matrix (see also jl^|). Thus, in the 
calculation of equations ^TM and (gJH, D xx = 1/6 = D ra , % = and D xy = 
(as Yq = 1 and J(y) — 0), so that Ao = 0, as expected for passive tracers, and 
D e = 1/6 + 6(8/945)Pe 2 . In the Bees & Croze theory D e and Pe are scaled with 
respect to the diffusivity scale V^/d r . Using the scaling (1 / 3)Vg / (2d r ) for random 
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diffusers in 3D, we recover D' e = 1 + (8/945)Pe /2 , the classical Taylor- Aris result 
for channels. 

The dispersion of passive tracers in turbulent channel flows has also been eluci- 
dated [§]. Elder derived the (dimensional) effective diffusivity K of passive tracers 
in turbulent open channel flow as K = 5.9 u T H, where u T is the friction velocity 
and H is the channel depth @; 0] . This open channel result applies equally to a 
closed channel with half- width H. Using the approximation U c /u T s» 5 log 10 Re [27j |. 
this result can be written as K ps 2.72vKe/ ln(Rc). To compare results across the 
laminar-turbulent transition (see later in figure I12j) . we use equations (|2.7|) and 
(|2.13|> to obtain a nondimensional turbulent diffusivity D e = K/Dq as a function 
of Pe, where Dq is the characteristic diffusivity scale. However, we stress that the 
diffusivity depends only on Re in the turbulent case: molecular diffusivity is negli- 
gible p . Notice also from the result of Elder how the effective turbulent diffusivity 
grows less sharply with Pe than in the laminar regime. Below, we will contrast the 
classical predictions for dispersion discussed in this section with our new results for 
biased swimmers. 



(b) Gyrotactic algae: a new, non-classical 'swimming dispersion' 

(i) Laminar downwelling flows: gyrotactic swimming strongly affects dispersion 

It is illuminating to consider how gyrotactic cells distribute across a channel in 
downwelling laminar flows, as this determines their streamwise dispersion. Figure 
2] displays the time-evolution of the cross-sectional cell concentration profiles in 
the wall- normal direction y for a selection of values of (Pe, 0). It is seen that an 
initially uniform concentration profile evolves to one focused at the channel cen- 
tre. This is what one would expect for swimming gyrotactic cells with orientations 
biased by combined gravitational and viscous torques 12; 13j]. At the population 



scale, the long-time concentration distribution is a result of the balance between 
cross-stream cell diffusion and biased swimming. In section (2£j), the cell concen- 
tration (normalised by its mean) is given by equation (gUP ): Y§ = c{y)/c = 
exp (/? J"p \q v /Dyv)ds) . (Recall subscripts m =G, F denote solutions of the GTD 
and FP microscopic models respectively.) The two models predict a qualitatively 
different dependence of concentration distribution on Pe and j3. Whilst the FP ap- 
proach allows for cells to diffuse more and focus less as the cells tumble in flows 
with large shear rates, the GTD model finds that with increasing local shear rate, 
a, the diffusivity of tumbling cells decreases faster than the decrease in cross-stream 
cell focusing. 

It has hitherto been unclear if GTD, mathematically more complicated than 
the FP approach, provides more accurate predictions. Thus we compare the DNS 
results and theoretical predictions from the two models. Figure [5] displays the long- 
time cell distributions from DNS contrasted with theoretical predictions using GTD 
and FP. It is clear that DNS profiles get broader with increasing /3/Pe, so profiles 
are shown for different values of this ratio (individual values of Pe and (3 are shown 
in the figure caption). GTD and DNS profiles agree very well for all the values 
of /3/Pe simulated. The FP model only agrees well with the DNS for large values 
of /3/Pe, where predictions coincide with those of GTD. The agreement is poor 
at low values (very poor for the sharply focused distribution f3/Pe= 0.002). It 
is possible to quantify the scaling of the breadth of the profiles with the ratio 
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Figure 4. Evolution of the DNS concentration profiles (not normalised) in the wall nor- 
mal direction y for gyrotactic cells in downwelfing channel flow with (Pe, /3) = (27, 1.34), 
(1675, 33.5), (670, 6.7), (1675, 3.35), from top to bottom as shown. At large times, the 
gyrotactic cells are observed to accumulate around the channel center. The profiles become 
more peaked as the ratio /3/Pe is decreased, see figure [5] 

/3/Pe. From asymptotic expressions for the ratio q v / D V q [l9j], it follows that GTD 
profiles can be conveniently approximated by Gaussian distributions of width j/o = 
[(2/A)(/3/Pe)] 5 (see appendix A). Clearly these predictions need also to be tested 
experimentally (see discussion), but the DNS results indicate that GTD is likely a 
superior model of gyrotactic response to flow than the FP approach. 
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Figure 5. Long-time concentrations profiles from DNS, normalised by the mean concentra- 
tion c, with (Pe, P, PfPe): (1675, 3.35, 0.002) (+), (670, 6.7, 0.01) (x), (1675, 33.5, 0.02) 
(A), (27, 6.7, 0.25) (□). Theoretical predictions from the GTD and FP models (lines, see 
legend), are compared with the simulations. GTD predictions agree very well with DNS 
profiles for all (3/Pe, but for FP agreement is poor for small values of this ratio. As in 
figure [4] DNS profiles broaden with increasing /3/Pe is increased. This is as predicted by 
GTD, where profiles can be approximated as Gaussians of width yo = [(2/A)(/3/Pe)] ' 5 , 
see text and appendix A. 

As for the passive case of figured we have quantified dispersion using the statis- 
tical measures: Ao(t), the drift above the mean flow; D e {t), the effective diffusivity; 
and 7(i), the skewness. These are plotted in figure H2 for Pe= 27, 670 and 1675 
at the fixed value of /3 = 33.5. In these simulations, statistically stationary values 
for Ao(t) and D e (t) are achieved for dimensionless times t ~ 1; steady values are 
reached earlier for larger Pe. In terms of the Bees & Croze dispersion theory, steady 
dispersion is achieved in the long-time limit when transient solutions to the mo- 
ment equations have died down: t 3> t\. The analysis of transient solutions with 
DNS will be carried out in a future study, but it is reasonable that gyrotaxis makes 
the approach to steady dispersion faster than for passive tracers. The skewness 
is negative and approaches zero (with 7 ~ £-0.49 f or p e _ 27; see supplementary 
materials), suggesting a distribution tending to normality. Bees & Croze predicted 
the power law decay 7 = 7o£~ 1 / 2 (the pre-factor 7o(/3,Pe) depends on gyrotactic 
swimming) with the same exponent as the passive case Q. 

The steady gyrotactic swimmer dispersion displays some very surprising fea- 
tures, evident in the data presented in figure^ For example, Ao, zero in the passive 
case, grows from a negative value to large positive values as Pe is increased. The 
effective diffusivity D e , on the other hand, shows a non-monotonic behaviour for 
increasing Pe (recall D e ~Pe 2 for passive tracers). We can qualitatively account 
for this behaviour considering the concentration distributions analysed above. Cells 
are biased to swim to the centre of the channel. Here only the torque due to gravity 
acts on cells, so cells swim upward at their mean swimming speed, which may be 
comparable with the mean flow speed for small Pe leading to Ao < 0. For large Pe, 
the upwards swimming speed is negligible. As cell accumulation at the centre of 
the channel increases with Pe due to gyrotaxis they drift more relative to the mean 
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Figure 6. Evolution of the dispersion of a population of gyrotactic swimming algae in a 
vertical down- welling laminar channel flow. The drift above the mean flow, Ao(t), effective 
diffusivity, D e (t), and skewness, 7(t) are displayed for Pe= 27, 670 and 1675, top to bottom 
as shown, for fixed /3 = 33.5. Uniquely, gyrotactic swimming algae have a non-zero drift 
above the mean flow, in contrast to the passive tracers in figure [3] Even more peculiarly, 
Ao < for Pe= 27. This is because of cell upswimming at mid-channel, see text. 

flow, and so Ao is an increasing function of Pe. Cell accumulation at the centre of 
the channel removes them from regions of high shear rate, eventually leading to a 
decrease in their effective axial diffusivity with increasing Pe. 

The Bees & Croze dispersion theory allows for the quantification of this intrigu- 
ing dispersion behaviour; the results are summarized in section (3st below. Prior 
to this, we shall consider time-dependent dispersion in turbulent flows. 

(ii) Turbulent downwelling flows: persistent but weaker swimming dispersion 

Here we describe the first DNS study of the dispersion of gyrotactic cells in 
turbulent channel flows. To compare turbulent and laminar results we first focus 
on downwelling flows. Simulations were performed for (Pe, Re)= (16570,2500), 
(28140, 4200) and (67000, 10000), with corresponding values of ft = 33.5, 55 and 67. 



10 2 

o 
< 

10 1 





Article submitted to Royal Society 



Dispersion of algae in channel flow 



15 



The long-time wall-normal cell concentration profiles for these turbulent flows are 
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Figure 7. Single realizations of the long-time concentration profiles (not normalised) in 
direction y for gyrotactic cells in downwelling turbulent flows for (Pe, Re)= (16570, 2500), 
(28140, 4200) and (67000, 10000), as shown. Gyrotaxis causes depletion of cells from regions 
close to the walls where shear is large, but the mean profile departs little from uniform 
concentration (expected of passive tracers). 

shown in figure [7J It is clear that the mean cell concentration is uniform (barring 
small fluctuations about the mean) for almost the entire channel width, except for 
regions close to the wall that are gyrotactically depleted of cells (gyrotactically de- 
pleted regions occupy only a small fraction (< 4%) of the channel width). Shear 
and advection experienced by a cell in turbulent flows are very different from the 
laminar case. The turbulent flow can be thought of as a time-averaged base profile 
on which are superposed turbulent fluctuations (the well-known Reynolds decom- 
position [27j ). The shear rate of the base profile is close to zero in the middle of the 
channel, and large at the walls, see figure [3J and deviates from the laminar case; 
this alone will lead to broader concentration profiles. On top of this, turbulent fluc- 
tuations perturb the flow, causing cell reorientation and advection. We can think 
of turbulence as endowing cells with an increased diffusivity (3o| acting to make 
gyrotactic accumulations in this downwelling case less pronounced. Only close to 
the walls is the impact of mean shear sufficient to cause significant cell depletion; 
an effect that increases with /3/Pe in a similar fashion to the laminar case. The 
small but measurable effects of gyrotaxis on the concentration profiles are reflected 
in the time-dependent dispersion measures for the turbulent case, plotted in figure 
[5] for the same values of (Pe, /?) considered above. We leave a detailed analysis 
of transients to a future study but note that, due to the increased diffusivity by 
turbulence, the approach to the limiting behaviour is faster than in the laminar 
case. Less obviously, the long-time dispersion retains a rich behaviour (observe the 
order of the curves in figure [7]). In particular, notice that cells have a nonzero drift 
Ao; this is due to local focusing of cells in downwelling regions of the fluctuating 
flow. The dispersion of gyrotactic swimmers is thus qualitatively distinct from that 
of passive tracers even in a turbulent flow. As Pe is increased to the maximum 
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Figure 8. Evolution of statistical measures (top to bottom, as for the laminar case) for gy- 
rotactic cells in a vertical downwelling turbulent channel flow for (Pe, Re)= (16570, 2500), 
(28140, 4200) and (67000, 10000), as shown. The long-time values of the drift above the 
mean flow, Ao and effective diffusivity, D e are plotted in figure 1121 Values of f3 are not 
shown for clarity, see table [T] 

value simulated, Ao decreases whilst D e increases, indicative of increased mixing 
by turbulence. 

(iii) Dispersion in upwelling flows 

The case of dispersion in flows directed vertically upwards (against the direction 
of gravity) is considered here briefly. The DNS results for the same values of Pe 
and j3 as the laminar downwelling case of section (SGJ are shown in Figure [§] The 
drift above the mean flow Ao is positive for small Pe, and grows more negative 
with increasing Pe. This behaviour is the result of the peculiar distribution of 
gyrotactic cells in the flow: cells in upwelling flow are biased to swim not to the 
channel centre, but to the walls [12]. Interestingly, accumulation and dispersion 
depend critically on the flow direction! The inset of figure H3 bottom, shows the 
normalised cell concentration c/c in the wall normal direction y, demonstrating the 
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Figure 9. Dispersion of gyrotactic algae in vertical upwelling laminar flows. The parameters 
are as for the downwelling laminar case (/3 = 33.5 and Pe = 27, 670, 1675). Top: the drift 
— Ao (inset: the positive drift for Pe = 27). Bottom: the effective diffusivity D e (inset: the 
normalised concentration profiles). 



wall accumulation, which becomes more peaked with decreasing /3/Pe (for more 
details, see supplementary materials). Thus, for fixed /3 and increasing Pe, cells 
focus at the walls and experience slower flow and less of the shear profile, leading 
to a decrease in both the drift, Ao, and the diffusivity, D e . The change of sign in 
Ao has a similar explanation as for the downwelling case. 

The upwelling turbulent case, presented in figure QUI was also investigated for 
the same parameter values as the downwelling case of section (^n]). We see that 
dispersion measures do not reach steady values: both Ao and D e grow monotonically 
with time for the duration of the simulation run. The inset to the top figure [TU] shows 
the cell concentration profiles displaying a thin band of gyrotactic accumulation at 
the walls. Cells that end up in this band are subject to strong dispersal by the large 
mean shear rate at the wall. Diffusion, whether due to turbulence or swimming, is 
not strong enough to balance the shear-induced migration towards the wall, as is 
the case for laminar flow, leading to non-steady dispersion over the course of the 
simulations. 

However, the upwelling dispersion results presented above are not realistic for 
swimming species that are negatively buoyant; accumulations of negatively buoyant 
cells at the walls can result in instability and the formation of downwelling flow and 
descending plumes near the walls [l2j]. These instabilities will differ in their extent 
in the laminar and turbulent upwelling cases, but we expect buoyancy driven wall 
flows to ensue after accumulation in both regimes. 
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Figure 10. Dispersion in upwelling turbulent flows for the same parameters as the down- 
welling case. Top: mean drift Ao (inset: cell concentration profile close to the wall). Bottom: 
effective diffusivity D e . Note that in the turbulent case neither of the statistical measures 
reaches a steady value. 



(iv) The effects of cell elongation 

So far we have assumed that the gyrotactic cells are spherical, setting the eccen- 
tricity parameter ao to zero. Here we present simulations obtained with a non-zero 
eccentricity, ao — 0.8, for the cases (Pe,/3) = (670, 6.7) (laminar) and (28140, 55) 
(turbulent). We chose the relatively large value of ao = 0.8 to emphasise the effect 
of eccentricity. Results for the effective diffusivity for downwelling flow are shown 
in figure QT] for the laminar and turbulent cases. The insets in the figure display 
the concentration profiles. We see that in the laminar case the effect of eccentricity 
is to broaden the profile, as observed in Bearon et al. 36]. This broader distribu- 
tion results in an increased value of effective diffusivity (cells sample more of the 
shear profile) . In the turbulent case of figure [TT] (bottom) there is a much smaller 
broadening effect. The effect of cell elongation on other statistical measures for the 
parameters considered here is marginal, see supplementary materials. If realistic 
values of biflagellate eccentricity (ao = 0-0.3 [1J|) are used, predictions for the 
dispersion of biased swimmers are not very different from those obtained here for 
spherical cells. 

(c) Long-time dispersion of gyrotactic cells as a function of Pe (Re) across the 

laminar-turbulent transition 

Having analysed the full time-dependence of gyrotactic cell dispersion, we con- 
centrate here on its long-time behaviour. In this limit, laminar DNS results can be 
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Figure 11. The effect of cell elongation on dispersion. DNS simulations of the effective 
diffusivity D e performed for spherical (ao = 0) and elongated (ao = 0.8) cells for down- 
welling flows. Top: laminar flow, with Pe= 670 and /3 = 6.7; bottom: turbulent flow, with 
Pe= 28140 and /3 = 55. Even for the large ao value considered, qualitative effects are not 
observed and significant quantitative differences are apparent only in the laminar case (see 
text and supplementary materials). 

compared with predictions from analytical theory for drift Ao and diffusivity D e 
as functions of Pe (for given f3) . The theory requires as inputs expressions for the 
mean swimming direction, q, and diffusivity tensor, D, obtained from microscopic 
stochastic models. We test two alternative microscopic models: Fokker- Planck (FP) 
and generalised Taylor dispersion (GTD). The models predict qualitatively differ- 
ent functional forms for the components of D as a function of the dimensionless 
shear a. Correspondingly, predictions for Ao(Pe) and D e (Pe) obtained with the two 
models also differ. 

The predictions from the analytical theory are shown as curves in figure [T^b . 
b, for fixed values of (3; the corresponding DNS results for selected Pe and the 
same values of fj are shown as points. It is clear that, for the GTD predictions, the 
agreement between theory and simulation for A is very good. The GTD prediction 
is that, for fixed /3, Ao increases with Pe, tending to a linear behaviour for large 
enough Pe. This is in contrast to passive tracers, for which Ao = 0. For small 
Pe the FP prediction coincides with that of GTD, but, for larger Pe, Ao tends 
to a constant instead of growing with Pe. The smaller the value of /3, the greater 
is the difference between the predictions of the two microscopic models (as for 
concentration profiles). A similar trend is seen for the effective diffusivity D ei shown 
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Figure 12. DNS results for long-time (a) drift, Ao, and (b) diffusivity, D e , as a function 
of Pe for gyrotactic cells in downwelling flows. Values of fi are: 1.34 (#), 6.7 (+), 33.5 
(□), 55 (a), 67 (♦). The DNS results are compared with predictions from the Bees & 
Croze dispersion theory, using the FP and GTD microscopic models, plotted as lines for 
the same j3 values as the DNS, see legend. Simulation results compare well with GTD 
predictions for all ft, but, for Pe> 1, they are incompatible with FP predictions for small 
p. Also plotted in (b) are the classical results for passive tracers. In laminar flows, Ao = 
and D e ~Pe 2 , in sharp contrast with our predictions for swimmers. In turbulent flow, we 
still have Ao > for swimmers, though it appears to decay to zero with increasing Pe. 
The turbulent D e for swimmers is close to the passive case; see text. 
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as a function of Pe in 112b. We see that for j3 = 33.5, both GTD and FP predict a 
diffusivity rising and then falling with Pe, consistent with DNS results. For smaller 
values of 0, however, the FP and GTD predictions are different: GTD predicts a 
decrease of D e at large Pe, while FP predicts a (second) rise. The DNS results are 
not compatible with this rise, but agree well with the GTD prediction. We remark 
that this dispersion behaviour is unique to swimming cells. It is due to gyrotactic 
bias and the ensuing accumulations that change the distribution of cells and thus 
their dispersion, as discussed above. Also shown in figure rT2"b is the diffusivity for 
passive tracers, which grows as D e ~ Pe 2 [5j, without the decrease at large Pe we 
predict for gyrotactic swimmers. 

In earlier sections we have seen how the effect of gyrotaxis on accumulation and 
dispersion is more subtle in turbulent flows. DNS results for long-time dispersion 
reflect these changes and are also summarised in figure [T^l An analytical theory of 
turbulent swimming dispersion has yet to be formulated, so we compare only with 
theory for passive tracers. It is seen for all simulations with (3 — 33.5, that the value 
of Ao, which was growing with Pe in the laminar regime, drops a little just beyond 
the transition to turbulence. On the other hand, the diffusivity D e , which was 
falling at large laminar Pe, suddenly acquires a sizeable value. The turbulent shear 
profile and fluctuations alter the distribution of cells. The effect of gyrotaxis is thus 
weakened but still measurable in the dispersion, which is similar but not identical to 
that of passive tracers. The fractional drift above the mean flow speed U = (2/3)U c , 
is given by S = lim t ^ tx) (l/C/)[dmi/(it — U] = (3/2)Ao/Pe; this is rather small for 
large turbulent Pe values. Nevertheless, it is remarkable that drift should not be 
zero in a turbulent flow, as it is for passive tracers. For very large values of Pe, 
we expect gyrotaxis to have practically no effect on dispersion: the time-averaged 
concentration profile will be well-approximated as uniform and Ao = 0. Indeed our 
results suggest that Ao — > for increasing Pe in the turbulent regime. 



4. Discussion 



We have studied the dispersion of gyrotactic swimming algae in channel flows by 
direct numerical simulations (DNS) and analytical theory. This is the first study 
to evaluate cell distributions and statistical measures of gyrotactic cell dispersion 
(drift above mean flow Ao, effective diffusivity D e and skewness 7) for flows on 
either side of the laminar-turbulence transition. We find unique cell accumulations 
and dispersion with non-zero drift and a non-monotonic diffusivity with Peclet 
number, Pe, with qualitative differences for upwelling and downwelling flows. The 
dispersion behaviour is remarkable and unique to gyrotactic swimming cells; passive 
tracers are transported at the mean flow rate (Ao = 0), the diffusivity increases 
with Pe (D e ^Pe 2 ) [5J, and dispersion is insensitive to channel orientation. In the 
laminar downwelling regime, simulation results were compared with the predictions 
derived from a recently formulated general analytical theory of swimming dispersion 
3 EH i applied here for the first time to channel geometry. The theory requires as 
inputs expressions for cell response to the local shear in a flow, determined from 
microscopic models. 

We have evaluated predictions based on two alternative models: the Fokker- 
Planck (FP) [H; Hi] and generalised Taylor dispersion (GTD) [IT ' 
proaches. Which model is more realistic has long been a matter of debate 
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we find that the DNS are in excellent agreement with analytical predictions based 
on GTD, providing the first evidence that it is much better than the FP model 
at describing swimmers in flows. As well as validating analytical predictions for 
laminar flow, the DNS allow us to explore the industrially relevant dispersion of 
gyrotactic algae in turbulent channel flows. In the turbulent regime the effects of 
gyrtoaxis (accumulations resulting in non-zero drift) persist, but are much more 
subtle. Effective diffusivity in downwelling turbulent flows is similar to that of pas- 
sive tracers. These are the first full direct numerical simulations of biased swimmers 
in turbulent channel flows; previous studies having concentrated on vortical flow 



37J or synthetic approximations of homogeneous isotropic turbulence [35|; |31| . 

The fact that swimming algae in channel flows distribute very differently to 
passive tracers has important practical consequences for the culture of swimming 
species. The most dramatic implications of our findings are for photobioreactors 
that operate using laminar flows. For example, in draft tube air-lift bioreactors, 
bubbles up a central draft tube (riser) mix and aerate cells, which then circulate 
down the channel formed between the draft tube and the surface of the reactor 
(downcomer; see figure 0J). The Bees & Croze [l[ analytical theory, confirmed by 
simulations, predicts that gyrotactic swimming algae will be focused more and 
more sharply at the center of the downcomer as Pe (the flow rate for fixed channel 
width, or Re) is increased. For example, considering a flow with Pe= 1675 (Re= 
250) and the swimming Peclet number (3 — 33.5 (H= 5cm), we predict that the 
concentration at the walls is a vanishingly small fraction of the mean, given by 
c/c w exp [— PeA/(2/3)] ~ 10~ 22 (using A = 2.2 for Chlamydomonas augustae] see 
approximations in Appendix A) . Non-swimming cells and molecular solutes, on the 
other hand, would be uniformly distributed across the tube, c/c = 1. As we can 
reasonably assume that the probability of cell adhesion to the walls is proportional 
to the concentration there, we predict that surface fouling by gyrotactic cells will 
be markedly reduced relative to non-swimming cells. Fouling can be a big problem 
in closed bioreactors because cells buildup can prevent light penetration and thus 
growth, and in extreme cases can clog reactor conduits |38j |. 

The peculiar dispersion of biased swimmers will also affect growth in biore- 
actors. Our results indicate that gyrotactic swimmers in a downwelling flow drift 
faster and diffuse less than passive tracers, which travel at the mean flow. We can 
make experimentally measurable predictions for the dispersion of Chlamydomonas 
augustae, the alga whose swimming parameters (reorientation time B, rotational 
diffusivity d ri swimming speed V s ) we have used in this study. For example, in 
the case with (3 = 6.7 and Pe= 670 (e.g. a realistic, small air-lift with H — 1 cm 
and U c = 1 cm/s) the fractional drift above the mean flow, 6 = (3/2)Ao/Pe, is 
estimated to be S — 0.48. In other words C. augustae cells drift about 50% faster 
than passive tracers, such as nutrients. For the nondimensional effective diffusivity 
we predict D e = 0.8: the axial diffusivity of cells is smaller than the cell diffusivity 
scale (V^/d r = 1.49 x 10 -3 cm 2 /s) in the absence of flow. To compare these pre- 
dictions with those for passive molecular solutes, we consider the dispersion of CO2 
(molecular diffusivity D t = 1.6 x 10~ 5 cm 2 /s), a vital nutrient for algae. In the 
flow under consideration, cells have Pe= 670, but CO2, with its smaller diffusivity, 
has Pe t =PeV*/(d r D t ) = 6.2 x 10 4 . Thus, the non-dimensional tracer diffusivity is 
D e ,t ~ (8/945)Pe 2 ~ 10 7 , from the Taylor-Aris result. Re-dimensionalising swim- 
ming and solute effective diffusivities by multiplying by diffusivity scales, we find 
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that the axial diffusivity of CO2 along a channel is « 10 5 times greater than that of 
C. augustae. This dramatic differential dispersion of cells and nutrients could have 
important consequences for swimming cell growth in reactors. 

In turbulent downwelling flows, which may also be realised in air-lifts, the ef- 
fects gyrotactic swimming are much less pronounced. Gyrotactic depletion is not 
as efficient as in laminar flows: the concentration at the walls is at best half of the 
mean concentration. We thus expect significant fouling in air-lifts under turbulent 
regimes. The influence of gyrotaxis on dispersion is also weaker than in the lam- 
inar case. For C. augustae in a turbulent channel flow with /? = 55 (H = 8.4cm) 
and Pe= 28140 (Re= 4200, U c = 5cm/s), the DNS predict a fractional drift of 
5 = 0.015. This is very small, and may be neglected for short channels. As shown 
in figure [T2"b. the effective diffusivity is of the same order as that of a passive tracer 
such as C02- 

The excellent agreement between DNS and predictions for dispersion obtained 
using generalised Taylor dispersion theory, provides a first confirmation that the 
swimming dispersion theory of Bees & Croze [l[ is robust. If they are to be useful in 
bioreactor engineering design, however, the theoretical and numerical predictions for 
dispersion need to be tested against experiments. Work is in progress to test GTD 
predictions for pipe flow with the alga Dunaliella 23] . It would also be interesting 
to test experimentally the numerical and analytical predictions for channel flow 
presented herein, exploring in particular the laminar-turbulent transition region. 
Since tubes are often arranged horizontally in bioreactor designs, it will be interest- 
ing to study the effect of tube orientation on biased swimmer dispersion. Croze et 
al. [39l | have carried out experiments on the dispersion of Chlamydomonas augus- 
tae in horizontal pipe flow for low Pe. They found a complex transport mediated 
by sheared bioconvection patterns and suggested that such cell-driven flows could 
alter the transition to turbulence. It would be interesting to study this transition 
experimentally and with the aid of simulations such as presented here, including 
the effects of cell buoyancy (especially for small Pe flows) [l9[ . 

More generally, our study could be extended to open channel flows, such as 
those present in pond bioreactors, channels, waterfalls and rivers. The swimming 
dispersion effects we have explored must also exist whenever biased swimmers are 
subject to shear in the ocean and lakes. Significantly, about 90% of species impli- 
cated in the formation of harmful algal blooms swim using flagella 40] . A recent 
study has shown that Heterosigma akashiwo is gyrotactic and could be responsible 
for thin layer formation 16]. However, the role of biased swimming in the dispersion 
and ecology of algae remains largely unexplored. 
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Table 2. Parameters used in the functional fits of solutions to the FP and GTD models 
used for our analytical theory predictions (from \nJ\j)- 



Appendix A. Microscopic model solutions 

(a) Numerical solutions for the Fokker-Planck and generalised Taylor dispersion 

models 

q(a) and D m (cr) predicted by the FP and GTD microscopic models have the 
following structure: 



q y(a) 





D, 



d%{?) Dyy(a) 










(Al) 



where we recall m =F represents solutions of the FP model and m =G those 
of the GTD model. Such solutions can be obtained using approximations via a 
Galerkin method |19j . It is convenient to fit such solutions with the following expres- 
sions: q x (a) = -P(a; a x ,b x ); qV(a) = -aP(a;a y ,b v ); D™(a) = P(cr; sl xx , h xx ); 
Dl>v{(j) = P{a;a.yy,byy);D x y(a) = -aP(a; a x «, b x ^). Here the rational function 
P(tr; a, b) is given by 



P(a;a,b) 



do + a 2 o- 



(A 2) 



and, for A = \/{2d r B) 



1 + b 2 o- 2 + b 4 a 4 ' 
2.2, the coefficients a and b are provided in table © below. 



(6) Approximate GTD profiles and width scaling 

Using asymptotic results derived in (l9j , we derive an approximation for the con- 
centration profiles predicted by the GTD model. These profiles are given by equa- 
tion (|2~^Tj) as c(y)/c = exp (/3 f^q* /D™)ds), where m = G for GTD predictions. 
For a <C 1 at leading order the GTD prediction asymptotes to q v (a) = — aJ\/\, 
where J\ is a known constant for A = 2.2. In the same limit, D^?(cr) = Ji/A 2 , 
so that (<7 y /£>^) (J «i = -o-X. For a > 1 at leading order, q«(cr) = -(2/3)A/ff 
and D'yy{a) = dx/a 2 , where di = 0.68 for A = 2.2. Thus q v / D V J w -aX is a 
reasonable approximation for all a. Recalling for channel flow a = (Pe//3 2 )y, we 
see that (3q v /Dq 1 « X(Pe//3)y. Inserting this expression in the equation above and 
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Box 
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128 x 129 x 128 
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10000 
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2tt x 2 x 4tt/3 


512 x 193 x 256 



Table 3. Large-scale and turbulent (or friction) Reynolds number for the simulations pre- 
sented here. The table also reports the simulation box size in units of H , the channel half 
width, and the resolution used for the simulations. 



integrating gives the Gaussian profile 

c(y) ^exp(y/y ) 2 , (A3) 

where 

y «[(2/A)(/3/Pe)] - 5 (A 4) 

is the width of the profile. This scaling is observed in the concentration profiles 
obtained from DNS; see main text. 



Supplementary materials 

SI. Numerical methods for the simulations 

The numerical code is an efficient pseudo-spectral solver for the three-dimensional 
incompressible Navier-Stokes equations (equations 2.3 or 2.6 in the main text) with 



a particle tracking algorithm [llj adapted to swimming micro-organisms. The ve- 
locity components of the fluid phase are expanded in both x (streamwise) and z 
(spanwise) direction with Fourier modes, whereas Chebyshev polynomials are em- 
ployed in the wall-normal y-direction. To advance the Navier-Stokes equations in 
time, we use a fourth order Runge-Kutta algorithm. Periodic boundary conditions 
are assumed in x and z, with no-slip at both walls y = ±1. More details about the 



code are given in [4l|. For the DNS with laminar flow, the flow is not simulated, 
but is taken to be the analytic parabolic Poiseuille solution to the Navier-Stokes 
equations . 

The micro-organisms are treated as point particles evolved by means of a La- 
grangian solver. The fluid velocities and the components of the velocity gradient 
tensor are interpolated from the Eulerian grid onto the particle positions using a 
tri-linear interpolation. The time advancement of the particle uses the same Runge- 
Kutta algorithm as the time-advancement of the fluid except for the noise which 
is integrated by an Euler-Marayuma method 42]. Equation (2.5) of the main text 
is solved using quaternions rather than angles |43J. We perform numerical simula- 
tions of turbulent plane Poiseuille flow at constant mass flux. The Reynolds number, 
domain size and resolution used for the results presented here are reported in ta- 
ble [3] Note that resolution and box size are chosen following the classic results for 



turbulent channel flow in 44: 45 



In the table, we report both the large-scale Reynolds number defined by Re — 
U c H/v where U c is the centerline streamwise velocity for the laminar flow of 
same mass flux, and the turbulent Reynolds number Re T — u T H jv. The latter 
is defined by the friction velocity u T = ^Jb\^ and the half-channel width, where 
cr w = v[dU / dy] wa \\ is the shear stress at the wall [27j |. 
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Note that the simulation domain is shorter than the distance travelled by the 
swimmers during dispersion; also, for large times, the area occupied by the cells 
becomes longer than our computational box. However, the computational domain 
can be considered long enough for the velocity correlation to be negligible at a 
distance of the order of the box length L x . Therefore, a very long domain for the 
swimmers can be assembled by means of copies of a single Eulerian computational 
box of length L x by using the method of repeated domains described in [l(| • 



S2. Evalutation of the analytical dispersion measures 

To evaluate the integrals in (2.20) and (2.22), it was simpler to numerically 
evaluate an equivalent system of steady differential equations (e.g. use the numerical 
solution to dY°/dy = fi(q y /Dyy)Y°, instead of (2.21)). The system of ODEs was 
solved using the Matlab bvp4c routine. Solutions to the system were found for using 
both the Fokkcr- Planck and GTD expressions given in Appendix A. For comparison 
with the analytical predictions, concentration profiles from DNS were normalised 
by dividing by the area under the profile (the number of cells) , evaluated using the 
Matlab trapz routine. 
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S3. Supplementary results 
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(b) 

Figure 13. Approach to zero of —7 (the negative of the skewness so we can use a log-log 
plot) for (a) passive tracers and (b) gyrotactic swimmers for the same values of Pe con- 
sidered in figures 3 and 6 of the main text, respectvely. We fit the passive case and the 
swimmer case with Pe= 27 with the power law 7 = "/ot~ a where 70 and a are constants. 
For tracers we find 70 = 0.19 and a — 0.51 and for swimmers 70 = 2.42 and a = 0.49. 
The pre factors are different, but both exponents are compatible, as predicted by Bees & 
Croze [l|] who calculated that the skewness for swimmers should approach zero as ~ t~ ' 5 
at long times, as it does for passive tracers Q. 
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Figure 14. Gyrotactic cell concentration profiles (not normalised) for upwelling laminar 
flows. Just as in the laminar case, gyrotactic accumulations become more pronounced as 
the ratio /3/Pe is decreased, but in the upwelling case these accumulations are at the walls. 
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Figure 15. Effect of cell elongation on dispersion measures for downwelling flows: (a) lam- 
inar (b) turbulent. Even for the unrealistically large value of eccentricity used («o = 0.8), 
the effect of elongation is small apart from the diffusivity in the laminar CcLSG, ELS discussed 
in the main text. 
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